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1 Introduction 

Unfortunately, most physical problems cannot be solved analytically, which means 
that we are not able to give their solutions in closed mathematical form. Thus, we are 
forced to adopt numerical methods to obtain results. But, how to keep these numerical 
solutions under control? In fact, one has to know — in some sense — what will emerge 
numerically, what will be the approximate result. Consequently, the first step, before 
starting the computer, will be to get a feeling for the output. First, one will analyze 
the dimension of the output; secondly, one will try to "guess" the order of magnitude 
of the desired output; and, thirdly, one will improve this first "guess." 

In these lectures we illustrate, on various examples, how to derive, even for rather 
complicated problems, quantities like upper bounds on energy eigenvalues analytically, 
and we improve our results by numerical methods. In both cases, we shall demonstrate 
how to use Mathematica, described as "a system for doing mathematics by computer" 
i- 

First, we will discuss how to handle dimensions and how to apply scaling arguments 
to the Schrodinger equation. After introducing the variational procedure, we calculate 
upper bounds on the energy eigenvalues of the Schrodinger Hamiltonian, using "paper 
and pencil" as well as Mathematica 3.0 [|l|. 

One of the main advantages of the use of Mathematica is the capability to represent 
results very easily in graphical form. This feature is illustrated in Fig. |I] for a particular 
D-wave eigenfunction of the hydrogen atom. 
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Figure 1: The eigenfu 
orbital angular mi 



2 Dimensional Reasonlin 

We all know that it is a horrible task t°^IA/liM|l^ 
it is customary to use a system of ''iiaiBKf'lwmti 



Using these units, we are able £ 
units of energy E. Glearly, 



Then 




Dim[£t] = Dihr 



nl = 0, relative 



ns. In particle ppysics, at least, 
led by 



uantity in, e.g., 



and, consequently, the dimension of time t is the~invel^e of the dimension of energy: 



Dim[t] 



1 



Bim[E] 



Bim[E 



-ii 



From the well-known relation between the mass m and the rest energy E of a particle, 

E = m c 2 , 
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the dimension of any mass m is 



Dim[m] = Dim[i?] 



Because of 

Dim[p] = Dim[mc] = Dim[m] , 
the dimension of any momentum p is also 



Dim[p] = Dim[i?] 



The dimension of any spatial coordinate x is obtained from 

Dim[pa;] = Dim[/i] = 1 



as 



Dim[x] 



Dim[p] 



Dim [E~ l ] 



We are now able to express any physical quantity in units of energy. For instance, let 
us calculate the dimension of the Schrodinger wave function ^ . The Schrodinger wave 
function \I/(x) in configuration-space representation is usually normalized according to 



J d 3 x^*(x)^(x) = 1 , 



which entails for the dimensions 



Dim[x 3 = 1 , 



hence, 



implying 



Dim[^ 2 



Dim [a; 3 ] 



Dim[tf] = Dim[£ 3 / 2 ] 



Let us apply, just as an example, this method of counting dimensions to the evaluation 
of the Fourier transform of 1/p 2 . As we know, the detailed calculation needs some care, 
and, involving a distribution, is not completely trivial. The quantity to be evaluated is 



dV 



ip-X 



Obviously, the only free parameter in the above integrand is x, and since the integrand 
contains the rotationally invariant scalar product p-x, the integral must be a function 
of merely the modulus |x| of the (spatial) coordinates x. Now, let us count dimensions: 



Dim 



= Dim[p] = Dim[£] = 



Dim[x] 
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Therefore, the integral has to be proportional to the inverse of |x|: 

/ d 3 p — — = — . 
J p 2 |x| 

The constant of proportionality A is now easily determined by applying the Laplacian 
A = V- V to both sides of this equation, and by using the rather well-known relations 

A t^t = -4tt5 (3) (x) 
|x] 

as well as 

J d 3 pe- ipx = (2tt) 3 5( 3 )(x) . 

This procedure results in 

A J d 3 p^P- = -J d 3 pe- ip * = -(2tt) 3 ^ 3 )(x) = A-^ = -4nA5 {3 \x) . 
From this relation we read off for the constant of proportionality A in the above ansatz: 

A = = 2tt 2 . 

47T 

Just by counting dimensions and by using general physical arguments, we find, without 
any difficult integrations in the complex plane, 

/d 3 p — = . 

J p x 

Along the lines sketched above, we will show how to solve rather complicated problems, 
just by applying general theorems and rules. 



3 The Schrodinger Equation 

Solving the Schrodinger equation is, in general, a somewhat tedious task. However, as 
long as we are merely interested in the dependence of the energy eigenvalues E on the 
parameters entering in the Schrodinger equation, we do not have to solve the equation 
at all. We can adopt physical arguments only, arguments like the scaling behaviour of 
the Schrodinger equation. 



3.1 Scaling Behaviour 

The time-independent Schrodinger equation for a generic power-law central potential, 
i.e., a potential of the form V(r) = a r n , which depends merely on the radial coordinate 
r = |x|, reads in configuration space, abbreviating the Laplacian again by A = V-V, 



■— + ar n \ tf(x) =EV(x) 
2 ii J 



(1) 
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with the reduced mass p of the two-particle bound system under consideration defined 
by 



m\ 7712 



77l\ + 7B2 

Let us start by scaling, just for fun, the coordinates x by some (arbitrarily chosen) scale 
parameter k, 

X = K p , 

which gives, because of 



(2) 



and 



K 



K n p n 



the scaled (time-independent) Schrodinger equation 

( A, 



+ an n p n \ V(kp) =E^(kp) 



\ 2/iK 2 j 

Multiplying now the above equation by 2 p k 2 leads to 



A p + 2paK z+n p n j * = . 

Apart from our pure amusement, what might be the deeper reason for this re-scaling? 
Obviously, we can choose k in such a way that the factor 2 pa disappears: 



„2+n 



1 



and thus 



2 pa 

\ l/(2+n) 



2 pa 



(3) 



For this particular choice of k, the scaled Schrodinger equation is cast into the form 

-A p + p n ) * = e * I , (4) 



with e some (dimensionless) number which has to be determined by solving the scaled 
Schrodinger equation. Note that this scaled Schrodinger equation is dimensionless. We 
have thus managed to separate physics from "pure" mathematics by the identification 

2 p k 2 E = 6 . 

Inserting Eq. (|3]) shows us how the energy eigenvalue E depends on the parameters a, 
p, and n entering in the Schrodinger equation: 



E = 


' a 2 ' 


l/(2+n) 

e 


[(2/i) n . 



(5) 



Now we are in the position to discuss the physical behaviour of the energy eigenvalues 
E for various power-law potentials without solving explicitly the Schrodinger equation. 
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Coulomb potential (n = —1): 

E = 2 fj, a 2 e . 

This energy eigenvalue is proportional to the mass and to the coupling constant 
squared. The energy eigenvalue E has to be proportional to the mass because the 
only parameter with dimension of energy entering in the Schrodinger equation is 
the mass /i. However, we cannot conclude from dimensional considerations only 
how the energy eigenvalue E depends on the (fine-structure) coupling constant a. 



Linear potential (n = 1): 



2 X 1/3 



In contrast to the Coulomb potential, this energy eigenvalue E is proportional to 
the inverse of the mass \x to the power 1/3. 

Logarithmic potential {n = 0): Because of the (easily to check) identity 

r n _ i 
In r = lim , 

n->0 n 

the logarithmic-potential energy eigenvalue E is obtained for n = in Eq. (|5]): 

E = ae . 



Note that in this case the energy eigenvalue E is independent of the mass. That 
means, that, for instance, the difference of energy eigenvalues of excited states is 
the same even for different masses \x. 

We have learned how to find the parameter dependence of the energy eigenvalues E of 
the Schrodinger equation without solving this differential equation explicitly. Thus we 
are able to check numerical solutions by calculating ratios which are independent of e, 
like for, e.g., the linear potential, 

Eh 

E-2 

In this way, we are able to get a feeling for the accuracy of any numerical calculation. 



Af2 



I/O 



3.2 The Variational Method: A Derivation of Upper Bounds 

Since we are interested merely in an explicit application of the variational method, we 
omit all mathematically subtle discussions and explanations here. A discussion of the 
history of inequalities and variational methods for eigenvalue problems can be found in 
Ref. 0, and some applications are given in Ref. ||. In Ref. J3j, a rigorous discussion 
of the problem can be found. For our purpose, in the frame of these lecture notes, it is 
sufficient to know that, for any given Hamiltonian H with eigenvalues Ek, k = 1, 2, . . ., 
upper bounds Ef. on Ek may be obtained by computing the matrix elements of H with 
respect to some suitably chosen set of orthogonalized trial functions which depend 
on some variational parameter A. These upper bounds can be improved by minimizing 
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the resulting upper bounds E k (A) on the "true" energy eigenvalues with respect to the 
variational parameter A. In order to obtain also upper bounds for the radially excited 
states, one has to diagonalize the corresponding energy matrix E^. The detailed steps 
of this general procedure are: 

1. Choose a set of mutually orthogonal trial states |\l/j(A)), i.e., (\l/j(A)|\l/j(A)) = 5ij. 

2. Determine the matrix elements of the Hamiltonian H using these trial functions: 

^.(A) = (^(A)|i?|^(A)} . 

3. Determine the roots of the characteristic equation 



det 



E ij (\)-E(X)6 lj =0. (7) 



4. These roots E(X) are upper bounds for any variational parameter A you choose. 

5. You may want to improve the upper bounds by minimizing -E'(A) with respect to 
A: 

dE(X) 



~». A r 



6. In this way, you obtain the minimal upper bound -Efc(A m in) on Ek (in your chosen 
sector of Hilbert space): 

Ek < E k (X min ) . (9) 
We are now in a position to apply the above concepts to various power-law potentials. 

3.2.1 Ground States 

Now, let us start by considering ground states. The prototype for all studies of this kind 
is the Coulomb problem. The Hamiltonian H under consideration here is simply given 
by 

H=f-°, (10) 
2 // r 

with the electromagnetic fine structure constant a. First, we have to choose some trial 
function. Let us take, for "didactic" reasons, the hydrogen ground-state eigenfunction 
(of course, you can also choose, e.g., Gaussian trial functions) 

#(r,A) = Ne~ Xr , A* = A > , 

with the variational parameter, A, and the normalization factor, N, to be determined. 
Secondly, calculate the normalization factor iV: 

d 3 x^*q = 1 = N 2 J d 3 xe- 2Xr = NHti Jdrr 2 e- 2Xr = NHn , 

where we have used 

11 




In Mathematica PJ, the procedure is the following: 
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MATHEMATICA 



Integration: 

In[l] := Integrate [r~n Exp [-lambda r] ,{r,0, Infinity}] 

Out[l]= If [Re [lambda] > && Re [n] > -1, lambda^ _1_n ^ Gamma [1+n] , 

/ — ; r dr ^ 

Jo e~ (lambda r) 

This output has to be interpreted in the following way: If 9ft A > and 9ft n > — 1, then 
the result of this integration reads A~ 1_n T(l + n), otherwise, the input is returned as 



Jo i^ dr ' 



which means that Mathematica is not able to find a solution. The normalization factor 
is thus 

A 3/2 



N = 



and the normalized trial function reads 

A 3 / 2 , 

In order to be on the safe side, we check the dimensions of this expression on both sides 
of this relation. As we have already discussed before, the left-hand side has dimension 

Dim[£ 3/2 ] . 

For the exponential exp(— Ar) to make sense at all, Ar has to be dimensionless, which 
requires 

DimfAr] = 1 , 

leading to 

D " n|A| = Dimfr] = ^ ' 

and, therefore, 

Dim[tt] = Dim[£ 3/2 ] = Dim[A 3/2 ] . 

Obviously, the dimensions, at least, are correct. Now we will show how the above steps 
are defined and calculated using Mathematica. We give the printout as it is returned. 



MATHEMATICA 



Defining the trial function: 

In[2] := psi [r_, lambda_] := Exp[-lambda r] 
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Integration: 

In[3] := 4 Pi Integrate [r~2 psi [r, lambda] ~2 , {r,0, Infinity}] 

Out [3]= 4tt If [Re [lambda] > 0, ~, /°V 2 lambda r r 2 dr] 

4 lambda 3 Jo 

This output has to be read as: If 3? A > 0, then the result is 

4A 3 ' 

otherwise, Mathematica returns the input, 



47T 



f 

Jo 



-2\r ^2 



r dr , 



indicating that Mathematica is not able to evaluate this input. The complete result is 

4tt „ /a 3 



N 



1 ~> N 



4 A 3 V 7T 

Next, we have to evaluate the expectation values of the kinetic energy with respect to 
our trial function. Since we are only interested in the energy of the ground state, which 
is always characterized by £ = 0, for our purpose the Laplacian A is effectively given by 



A= -f- 



dr 2 r dr 



Its expectation value thus reads 

J d 3 xty*(r,\) Atf(r,A) 



7r J \ dr z r dr ) 

A 3 7 

— Air J dr (r 2 A 2 - 2 r A) e 

r(2)" 



-2 Ar 



A2 r(3) 

. (2 A) 3 



4A 



-2A 



(2 A)- 



= -A 2 . 

Again we will check the dimensions and we will omit from now on the symbol Dim[. . .]: 

x 3 ^\y ~E~ 3 E 3 / 2 E 2 E 3 / 2 = E 2 ~\ 2 . 
x 2 

Fortunately, all dimensions are O.K. The expectation value of the kinetic energy is thus 



(12) 



Already knowing how to interpret Mathematica output, below we just give the output. 
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MATHEMATICA 



Defining the trial function: 

In[4]:= psi [r_, lambdaj := lambda" (3/2) /Sqrt [Pi] Exp[-lambda r] 

Effective Laplacian for vanishing orbital angular momentum: 

In[5]:= g[r_, lambdaj := D [psi [r, lambda] ,{r, 2}] +2/r 
D [psi [r, lambda] ,{r,l}] 

Integration: 

In [6] := 4 Pi Integrate [r~2 psi [r, lambda] g [r, lambda] ,{r,0, Infinity}] 
lambda 2 \ 



Out [6]= 47T 



47T 



The expression D [psi [r , lambda] , {r ,n}] is the nth derivative of \l/(r, A) with respect 
to r. The expectation value of a power-law potential V(r) = ar n reads 

\3 °° 

J d 3 x^*(r,X)V(r)^(r,X) = — An a J drr n+2 e~ 2Xr 

A3 T{n + 3) 
vr (2A)™+ 3 ' 

that is, 



(V(r)) = I d 3 x **(r, A) V(r) #(r, A) = 4 a A 3 T ^ + 3) 



(2 A) 



n+3 



(13) 



For checking again the dimensions, we need the dimension of the coupling constant a: 

V = ar n ~ aE- n ~ £ 

implies 

a ~ £/ 

Without surprise, we realize that the dimensions in Eq. ( |13D are O.K.: 

a\- n ~ E n+1 E- n = E . 
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MATHEMATICA 



Integration: 

In[7]:= Integrate [r~(n+2) Exp[-2 lambda r] ,{r,0, Infinity}] 



Out [7]= 2 3 n lambda 3 11 Gamma [3+n] 



Combining Eqs. (|12[ ) and (|13|) yields the variational energy 

A 2 a i> + 3) 
E(A) "2^ + 2l2AF- 



(14) 



This energy is always an upper bound to the "true" energy for any A > 0: E trnc < E(X). 
We are easily able to improve this upper bound -E"(A) by determining that value A m i n of 
the variational parameter A which minimizes E(X); clearly, A m i n is then obtained from 



dE(X) 
dX 



. 



The derivative of Eq. (|1J) with respect to A entails the requirement 

dE(X) X I> + 3) 



dX 



an ■ 



which is solved by 



A, 



(2A)™ +1 

/ . Q N-|l/(n+2) 

anjjii {n + 6) 







yn+i 



(15) 



which value, when re-inserted into Eq. fli~4]), entails, in turn, the improved upper bound 
E var = E(X min ) = -(-) ±— -L 1 + - . (16) 



n 



MATHEMATICA 

Defining the function E(X): 

In [8]:= e[lambda_] := lambda~2/(2 mu) + a/2 Gamma [n+3] /(2 lambda) "n 

Differentiation with respect to the parameter X (note that D [e [lambda] , {lambda, 1}] 
is equivalent to D[e [lambda] , lambda] ): 

In [9] : = D [e [lambda] , lambda] 

_ _ lambda -i-n -1-n n 

Out [9]= - 2 1 n a lambda 1 n n Gamma [3+n] 

mu 
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Solving the equation for A min : 

r lambda -1-n -1-n ,- 

In [10]:= Solve [ - 2 1 n a lambda 1 n n Gamma [3+n] == 0, lambda] 

mu 

Out [10]= {lambdas (2 _1_n a mu n Gamma [3+n] ) 1 



Calculation of E (A m i n ): 

_._ 1 
In[ll]:= e[(2 1 n a mu n Gamma [3+n] ) 2+n ] 

r n (2 a mu n Gamma [3+n] ) 2+n 

Out [11]= + 2 1 11 a 

2 mu 

/ 1 \ _n 

Gamma [3+n] ( 2 1 11 a mu n Gamma [3+n] ) 2+n 



Note: The command PowerExpand [expr] expands all powers of products and powers. 
% stands for the last expression (here Out [6] ). 

In [12]:= PowerExpand [%] 

-4-3n 2 _ n _ n 2 
Out [12]= 2 2+n a 2+n mu 2+n n 2+n (2+n) Gamma [3+n] 2+n 

Out [12] is the same result as in Eq. fll6|) . Let us apply Eq. (|16|) to typical potentials: 
Coulomb potential: Inserting here a = —a and n = —1, we obtain the upper bound 



a 2 \x 



y true 



Interestingly, the expression on the right-hand side is precisely the ground-state 
energy of a hydrogen atom, in other words, here the sign of equality applies. The 
obvious reason for this clearly is that we have used the hydrogen wave function of 
the ground state as our trial function. For a = 1 and fi = 1, the numerical value 
for E is E = —0.5. It is very simple to find the minimum of Eq. (|14D numerically 
with the help of Mathematica. 



MATHEMATICA 
Defining the function to be minimized: 

In[13] := e[lambda_,n_,a_,mu_] := lambda~2/(2 mu) + a/2 Gamma[n+3] /(2 
lambda) ~n 

Finding the minimum (the starting point for the variational parameter A is A = 0.5): 
In [14] := FindMinimum [e [lambda, -1,-1,1] , {lambda, . 5}] 
Out [14]= {-0.5,{lambda^l.}} 



This has to be understood as E miri = 



—0.5 at the point A min = 1. 
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Linear potential: V(r) = ar. Inserting n = 1 into Eq. ([16]) entails 




Let us rewrite this result in the form of Eq. (|6]) : 

= 24764 y • (i7 > 

The true ground-state energy E true is given by the first zero of the Airy function 

!l, 

/a 2 \ 1/3 
-E'trnp = 2.3381 — , 
\2fj,J 

which implies a surprisingly small relative error of our crude upper bound E var : 

6%. 

-'-'true 

We have calculated, without solving the Schrodinger equation, the ground-state 
energy for the linear potential to unexpectedly good approximation analytically. 



E VSLT -^true 



MATHEMATICA 
Defining the function to be minimized: 

In [15] := e[lambda_,n_,a_,mu_] := lambda~2/(2 mu) + a/2 Gamma[n+3] /(2 
lambda) ~n 

Finding the minimum (the starting point for the variational parameter A is A = 0.5): 
In [16] := FindMinimum [e [lambda, 1,1,1] , {lambda, . 5}] 
Out [16] = {1 . 96556 , {lambda^l . 14471}} 

These results may be verified by inserting n = a = fi = 1 into Eqs. ( |TTD and ([IS]) . 
Sometimes it is illustrative to compare plots. Letting 2 fi = 1, we shall plot, first, 
i?true = 2.3381 a 2//3 , then, E v3iV = 2.4764 a 2//3 , and, finally, we shall combine these 
two plots into one single plot with the help of the Mathematica command Show. 



MATHEMATICA 
Defining the functions to be plotted: 
In [17]:= etrue [a_] := 2.3381 a" (2/3) 
In [18]:= eupper[a_] := 2.4764 a" (2/3) 
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Plotting E truc (plot 1 J: 

In [19]:= plotl = Plot [etrue [a] ,{a,0,3},AxesLabel->{"a\ "E"}, 
TextStyle->{FontSlant->" Italic" ,FontSize->14}] 

AxesLabel defines what shall be written on the axes; TextStyle defines the font type 
(here italic) and the font size (here 14 pt). 



Figure 2: The 
4 
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However: how is 
drawn the plot of 
a ps file. For this 



thi 



true energy E truc for the ground state ofthe^hnear potential 



s figure — and all the others — imported into Tj^K? After you have 
function using the command above, you must save the graphic as 
pprpose, youwrite in Mathematica, after the appearence of the plot: 



the. 



Display ["etrue .ps Vplotl , "epsi" , ImageResolution->300] 



etrue . ps is the name yoi^ give to tlje output jfiie, plot 1 is^he^ame ^f tSe graphics you 
gave to the plot (see above") and epsi is the necessary format for the graphics output. 
The command ImageResolution can be omitted, it depends on the resolution for the 
graphics you prefer. As the next step, you include your ps file in TpX with the help of 
the following commands: 



\begin{f igure} [h] 
\begin{center} 

\psf ig{f igure=etrue.ps, scale=l} 

\caption{etrue} 

\end{center} 

\end{f igure} 
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The command scale enables you to resize the picture. For instance, scale=0 . 5 would 
make this figure precisely half as large as the original one. The corresponding style file 
is included by the TpX command \documentstyle [epsf ig, . . . ] . In order to view the 
graphics, the dvi file has to be converted into a ps file. 
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Plotting E vlir (plot2j: 

In [20]:= plot2 = Plot [eupper [a] , {a, , 3} , AxesLabel->{"a" , "E"} , 
TextStyle->{FontSlant->" Italic" ,FontSize->14}] 



E 



4 



Figure 3: Variational upper bound E v3iI on the energy for the^ 
potential 



3 

Combining plotl 
In [21] := Showfpl 




nd state of the linear 



and plot2: 
otl,plofe2] 



1 

3.2.2 Radial Eieftations 



Up to now we have discussed, ground- statje energies only, ^Nfegt, we would like to show 
how to determine upper bounds on energies of radially excited states. A more detailed 
explanation of this procedure may be found in, e.g., Ref. [||. We illustrate this method 
by determining the ground-state energy and the energy of the first radial excitation of 
a system bound by some linear potential. First, we have to choose some trial functions. 
For the ground state, we take again the hydrogen eigenfunction 



*o(A,r) 



A 3/2 



-Xr 



A* = A > 



and, for the first radially excited state, we use 

A 3/2 



*i(A,r) 



/37T 



(3-2Ar)e 



-Xr 



A > . 



;is) 



Since the radially excited states are states with nodes in the wave function, we have to 
construct trial functions with this feature. The first radially excited state has one node, 
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thus the trial function has to have one zero, as is the case in Eq. (]XSf) . This trial function 
of the first radially excited state is constructed to be orthogonal to the eigenfunction of 
the ground state: 

(* |*i) = 0. 

To begin with, it's straightforward to check the orthogonality employing Mathematica. 

MATHEMATICA 
Defining the ground-state trial function \l/o-' 

In [22] := psiO [lambda_ , r_] : = Sqrt [lambda~3/Pi] Exp [-lambda r] 

Defining the trial function for the hrst radially excited state: 

In[23]:= psil [lambda_,r_] := Sqrt [lambda~3/ (3 Pi)] (3-2 lambda r) 
Exp [-lambda r] 

Checking the normalization of ^q: 

In[24]:= 4 Pi Integrate [r~2 psiO [lambda, r] "2, {r,0, Infinity}] 

r n 1 

Out [24] := 4 Pi 

4 Pi 

Checking the normalization of ^i: 

In[25]:= 4 Pi Integrate [r~2 psil [lambda, r] "2 , {r,0, Infinity}] 

r h 1 

Out [25] := 4 Pi 

4 Pi 

Checking the orthogonality of the wave functions \&o and ^i: 

In[26]:= 4 Pi Integrate [r~2 psiO [lambda, r] 
psil [lambda, r] ,{r,0, Infinity}] 

Out [26]= 4 Pi 



18 



The expectation value V 00 = (ty \ ar \ ^f ) of the linear potential: 

In[27]:= 4 Pi Integrate [r~2 psiO [lambda, r] a r 
psiO [lambda, r] , {r , , Infinity}] 

3 a 

Out [27]= 4 Pi 



8 lambda Pi 

The matrix elements Vbi = Vio = (^ol ar l^i) of the linear potential: 

In [28]:= 4 Pi Integrate [r~2 psiO [lambda, r] a r 
psil [lambda, r] ,{r,0, Infinity}] 

r -. , Sqrt[3] a N 

Out [28]= 4 Pi (-— ) 

8 lambda Pi 

The expectation value V n = ar \^fi) of the linear potential: 

In [29] := 4 Pi Integrate [r"2 psil [lambda, r] a r 
psil [lambda, r] ,{r,0, Infinity}] 

Out [29]= 4 Pi 



8 lambda Pi 

Defining the Laplacian of the ground-state trial function ^ : 

In[30]:= laplacepsiO[lambda_,r_] := D [psiO [lambda, r] ,{r,2}]+2/r 
D[psiO [lambda, r] ,r] 

Defining the Laplacian of the trial function for the first radially excited state: 

In [31] := laplacepsil [lambda_,r_] := D [psil [lambda, r] , {r, 2}] +2/r 
D [psil [lambda, r] ,r] 



The expectation value T 00 = ( \I/ C 



P 2 



2/i 



/ of the kinetic energy: 



In [32] := 4 Pi Integrate [r"2 psiO [lambda, r] (-laplacepsiO [lambda, r] / (2 
mu) ) , {r , , Infinity}] 

_ _ lambda 2 

Out [32]= 4 Pi 

8 mu Pi 

,2 



The matrix elements T m = T w — ( ^ c 



P 



2/i 



) of the kinetic energy: 



In [33] := 4 Pi Integrate [r"2 psiO [lambda, r] (-laplacepsil [lambda, r] / (2 
mu) ) , {r , , Infinity}] 

r n lambda 2 
Out [33]= 4 Pi 



4 Sqrt[3] mu Pi 
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The expectation value T u = ( 



P 2 



2/i 



) of the kinetic energy: 



In[34] := 4 Pi Integrate [r~2 psil [lambda, r] (-laplacepsil [lambda, r] / (2 
mu) ) , {r , , Infinity}] 

n _, 7 lambda 2 

Out [34]= 4 Pi 

24 mu Pi 

With the matrix elements calculated above, we are able to solve Eq. (0). However, it is 
much simpler to determine the energy eigenvalues using Mathematica. We shall show 
this analytically as well as numerically. 

MATHEMATICA 

The energy matrix elements: 

lambda 2 3 a 

In[35] := eOO [lambda_,mu_,a_] := + 



2 mu 2 lambda 

7 lambda 2 5 a 

In[36] := ell [lambda_,mu_,a_] := + 



In[37] := elO [lambda_,mu_,a_] : = 



6 mu 2 lambda 

lambda 2 a 



•y/3 mu 2 lambda 
Defining the energy matrix Eij(X) as a function of A, fi, and a: 

In [38] := ematrix [lambda_,mu_, a_] : = 
{{eOO [lambda, mu, a] , elO [lambda, mu, a] } , 
{elO [lambda, mu, a] ,ell [lambda, mu, a] }} 

Defining the energy eigenvalues E(\) as a function of A, \x, and a: 

In[39] := eeigen[lambda_,mu_,a_] := Eigenvalues [ematrix [lambda, mu, a] ] 

Analytic evaluation of the eigenvalues: 

In [40] := eeigen [lambda, mu, a] 

Out [40]= {(5*lambda~4*mu + 12*a*lambda*mu~2 - 
2*lambda*mu*Sqrt [4*lambda~6 - 6*a*lambda~3*mu + 

9*a~2*mu~2] )/(6*lambda~2*mu"2) , (5*lambda~4*mu + 12*a*lambda*mu~2 + 
2*lambda*mu*Sqrt [4*lambda~6 - 6*a*lambda~3*mu + 
9*a~2*mu~2] ) / (6*lambda~2*mu~2) } 

Simplifying the above expression: 

In [41]:= FullSimplify[ /.] 

Out [41]= {(5*lambda~3 + 12*a*mu - 2*Sqrt [4*lambda"6 - 6*a*lambda~3*mu 

+ 9*a~2*mu~2] ) / (6*lambda*mu) , (5*lambda~3 + 2*(6*a*mu + 

Sqrt [4*lambda~6 - 6*a*lambda~3*mu + 9*a~2*mu~2] ))/(6*lambda*mu)} 
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Numerical evaluation of the eigenvalues E for A = 1 GeV, \i = 0.5 GeV, a = 1 GeV 2 : 
In [42] := eeigen [1 , 1/2 , 1] 

Out [42]= {(11 - Sqrt[13])/3, (11 + Sqrt[13])/3} 

Numerical value (the Mathematica command N ( expr) yields the numerical value of the 
expression expr): 

In [43] := N[°/„] 

Out [43]= {2.46482, 4.86852} 

The result tells us: the variational upper bounds for the energy eigenvalues are given by 

^uppor(lS) = 2.46482 GeV , 

Supper (2S) = 4.86852 GeV , 

which have to be compared with the true values (determined by the first two zeros of 
the Airy function) 

£ truc (lS) = 2.33811 GeV , 
£ truc (2S) = 4.08795 GeV . 

For the ground state, labelled IS in usual spectroscopic notation, the relative error of 
our variational upper bound -& up per is 

-^uppcr(lS) — gtruejlS) _ g ^ y 
-Etrue(lS) 

which is, obviously, slightly better than the previous relative error of 6 %. Accordingly, 
increasing the matrix size — here from 1 x 1 to 2 x 2 — improves our upper bound. For 
the first radially excited state, labelled 2S in usual spectroscopic notation, the relative 
error of our variational upper bound -Supper is 

-Supper (2S) — -^truc(2S) _ ^ ^ 
Strue(2S) 

which is, admittedly, a little bit large. However, we are able to improve these results by 
applying Eqs. @ and In this way, we find the minimum of the energy eigenvalues 
covered by our chosen set of trial functions. 

However, varying the variational parameter A usually destroys the orthogonality of 
trial functions attributed to different radial excitations. The characteristic equation is, 
in general, not given by Eq. (^). (A detailed discussion of these troubles can be found in 
Ref. ||.) In view of this, we focus our attention in the following to the ground state IS, 
where this problem does not show up. 

Applying the minimization procedure described at lengthy at the beginning of this 
section, we try to optimize our previously obtained upper bounds on the ground-state 
energy of the linear potential by minimizing the corresponding eigenvalue E(X) of our 
energy matrix Eij. This eigenvalue is extracted from the above analytical results using 
the Mathematica command Part [expr, i] which returns the ith part of the expression 
expr. 
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MATHEMATICA 



Defining the ground-state eigenvalue as function of A for u, = 0.5 GeV and a = 1 GeV 2 : 

In [44]:= eOOeigen [lambda.] : = 

Part [Eigenvalues [ematrix [lambda, 1/2 , 1] ] ,1] 

Calculating the ground-state energy eigenvalue analytically: 

In [45]:= eOOeigen [lambda] 

Out [45]= (6*lambda + 5*lambda~4 - lambda*Sqrt [9 - 12*lambda~3 + 
16*lambda~6] ) / (3*lambda~2) 

Finding the minimum for the lS-state energy (the starting point for A is A = 0.5 GeV): 
In [46] := FindMinimum[°/ , {lambda, 0.5}] 

We might expect that the correct minimum may be found employing the Mathematica 
command FindMinimum [eOOeigen [lambda] , {lambda , . 5}] . Unfortunately, for very 
mysterious reasons, Mathematica returns a wrong result. Because of this, we carefully 
split the calculation into parts. Presumably, Mathematica is not able to simultaneously 
calculate eigenvalues, extract parts of a matrix, and find minima of some expression. 

Out [46]= {2.4322, {lambda->0 . 665633}} 



The new result, improved by minimization with respect to the variational parameter A 
within the set of upper bounds obtained before, thus reads 

£ var (lS) = 2.43220 GeV 

at the point 

A min = 0.665633 GeV . 

Evidently, we succeeded in reducing the corresponding relative error significantly, viz., 
to 

-5>ar(lS) - gtrue(lS) = 3 g % 
-Etrue(lS) 

As mentioned before, minimization with respect to the variational parameter leads, 
in general, to different values for this variational parameter in the ground state and in 
the radially excited states: Aj 7^ Xj. Therefore, in general, all these states are no longer 
necessarily orthogonal: 

(v| M A,) vl/ ; ( A; )) / % . 

In this case, the characteristic equation for our eigenvalue problem generalizing Eq. (|7]) 
reads 



det 



. 



_<*,(A,) \H\y j (\ J ))-E(y i (\ t )\y j (\ J )}_ 

In the next subsection we will show how to increase the matrix size — and therefore 
the accuracy of the computed upper bounds on the Schrodinger energy levels — further. 
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3.2.3 Laguerre Bounds 

As mentioned before, the crucial step in obtaining accurate upper bounds is the choice 
of convenient trial functions. Here we shall work in a basis which involves the so-called 
generalized Laguerre polynomials U}J [|7| , specific orthogonal polynomials, defined by 
the power series || 

4 7) (*) = E(-i) J 

j=0 

and normalized according to pf 




(7) M r(7 + fc + i) 

X) — Tj Okk' 



Consequently, introducing now even two variational parameters, namely, one, A, with 
the dimension of mass as well as a dimensionless one, /3, our choice for the trial function 
, 0fc,im( x ) with orbital angular momentum I and projection m, reads in coordinate space 

^Mm(x) =AT|x|^- 1 exp(-A|x|)4 7) (2A|x|)^ m (fi) , (19) 

where normalizability constrains the variational parameter A to positive values: A > 0. 
Here, 3^m(^) label the spherical harmonics for angular momentum i and projection m 
depending on the solid angle Q; by convention, they are orthonormalized according to 

an y\ m (o) y Vm , (n) = s u < <W . (20) 

The proper orthonormalization of our ansatz (|i~9l) fixes the parameter 7 necessarily to 
the value 7 = 2 £ + 2 (3 and determines the normalization constant TV: 



2A)2Wlfc! exp(-A |x|) 4^(2 A |x|) ^(fl) 



^ r(2£ + 2/5 + A; + l) 1 1 
satisfies the normalization condition 

d 3 X-0^ m (x) i>k',l'm'{x) = 8kk' hi' $mm' 



Rather obviously, normalizability constrains the second variational parameter, /3, too, 
namely, to a range characterized by 2 (3 > —1, i.e., to the range (3 > — \. For simplicity, 
we set our mass scale by choosing m\ = m 2 = 1 GeV and fix the variational parameters 
to the values A = 1 GeV and (3 = 1. We illustrate the general procedure by discussing, 
for comparison, again the linear potential V = r, for simplicity with slope a = 1 GeV 2 . 
At least up to a 4 x 4 matrix, all calculations may be performed analytically; we leave 
this to the reader as some exercise, and follow our procedure by applying Mathematica: 

• define the chosen set of trial functions ^^(x); 

• compute the matrix elements of the Laplacian, in other words, the kinetic energy; 

• compute the matrix elements of the radial coordinate r, i.e., the potential energy; 

• determine the eigenvalues E(\) of the resulting matrix i%(A) of the total energy; 

• compare the eigenvalues for different matrix sizes in order to observe, hopefully, 
convergence. 
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MATHEMATICA 
Defining the trial function ipk/mfa): 

In[47]:= psix [k_, l_,m_, r_] := Sqrt[2~(2 1+3) k ! /Gamma [2 1+3+k]] r~l 
Exp [-r] *LaguerreL [k, 2 1+2,2 r] *SphericalHarmonicY [1 ,m, theta,phi] 

Since we discuss S waves only, we define the trial function for £ — m — 0: 

In[48]:= psi[k_,r_] := psix[k,0,0,r] 

Defining the Laplacian A^(r) operating on states with i = (S waves): 
In[49]:= delta[k_,r_] := D[psi [k,r] ,{r,2}]+2/r D[psi [k,r] ,{r , 1}] 
Defining the integrand ips{r) Aipk(r): 
In[50] := intks [k_,s_,r_] := psi[s,r] delta[k,r] 

oo 

The matrix elements j drr^ s {r) (-A^(r)) ofthekinetic -energy operator T = — A.fl 
o 

In [51]:= kinen[k_,s_] := -4 Pi Integrate [r" 2 
intks [k , s , r] , {r , , Inf inity}] 

oo 

The matrix elements J drr 2 ip s {r) ripk{r) of the potential-energy operator V(r) = r: 

o 

In [52] := poten [k_, s_] := 4 Pi Integrate [r "3 psi[s,r] 
psi[k,r] ,{r,0, Infinity}] 

The matrix elements of the total energy: 

In [53] := toten[k_,s_] := kinen [k, s] +poten [k, s] 

Here we construct the matrix built by the above matrix elements, using the command 
Table. Since counting of the matrix indices starts at 0, we redefine this matrix in order 
to obtain for x = 1 a 1 x 1 matrix, and so on. 

In [54]:= totenmat [x_] := Table [toten [k, s] , {k,0,x-l} , {s ,0 ,x-l}] 

Defining now the function eeigen [x] with the help of the command Eigenvalues [M] 
which gives the eigenvalues ofanxxx matrix M, in other words, diagonalizes any xxx 
matrix (for our purposes, it has to be applied, of course, to the matrix totenmat [x] of 
the total energy): 

In [55] := eeigen [x_] := Eigenvalues [totenmat [x] ] 

1 Since we have chosen both particle masses to be mi = ni2 = 1 GeV, twice the reduced mass /x is 
also 1 GeV. 
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Eigenvalue of the lxl energy matrix: 
In [56] := eeigen[l] 

Out [56]= jj} 

Eigenvalues of the 2x2 energy matrix: 
In [57] := eeigen[2] 



Using the command N [°/„] for numerical evaluation of the last output: 
In [58] := N[°/„] 

Ground-state energy and energy of the first radial excitation numerically: 
Out [58]= {2.46482,4.86852} 
Eigenvalues of the 3x3 energy matrix: 
In [59] := eeigen[3] 



Using the command N [°/J for numerical evaluation of the last output: 



Ground-state energy and energies of first and second radial excitations numerically: 
Out [60]= {4.5,2.35425,7.64575} 

In general, the eigenvalues of a 5x5 matrix cannot be found analytically; we determine 
them numerically: 

In [61] := N[eeigen[5]] 

Ground-state energy and energies of the first four radial excitations: 
Out [61] = {2 . 34136 , 4 . 13334 , 5 . 72535 , 8 . 11424 , 15 . 519} 

The computation of the eigenvalues of the 10 x 10 energy matrix consumes already a 
noticeable amount of computer time: 

In [62] := N[eeigen[10]] 

Ground-state energy and energies of the first nine radial excitations: 



{2 . 33812 , 4 . 08858 , 5 . 53209 , 6 . 83859 , 8 . 14892 , 9 .91409 , 12 . 195 , 14 . 096 , 17 . 146 , 




Out [59]= {-, 5 - V7, 5 + V7 



In [60] := N[°/„] 



Out [62] = 



49.7026} 
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Table [T] compares the dependence of the accuracy of the obtained bounds on the energy 
levels of ground state and first radially excited state on the size of the energy matrix. 

Table 1: Comparison of relative errors of the energy eigenvalues for varying matrix sizes 



Matrix Size 


IS State 


2S State 


1 x 1 


6% 




2x2 


5 % 


19 % 


3x3 


0.7% 


10 % 


5x5 


0.1 % 


1 % 


10 x 10 


4 x 10" 4 % 


2 x 10- 2 % 



Obviously, the accuracy obtained is very impressive. One should keep in mind that the 
above procedure is applicable not only to potentials V(r) different from the power-law 
form V = r n but also to differential operators different from the nonrelativistic p 2 /2 m. 
In this way one is able to determine even the energy eigenvalues of more sophisticated 
Hamiltonians like the semirelativistic one, which incorporates the square-root operator 
y/p' 2 + m 2 of the relativistic kinetic energy J7|. Moreover, it might be of interest to note 
that up to and including 4x4 matrices, the diagonalization of the energy matrix E can 
be performed analytically. This fact enables one to gain control over purely numerically 
obtained results. However, one should never forget that the numerical values computed 
here represent always only upper bounds on the true energy eigenvalues: E tTue < E(X). 

In summary, we have tried to show how (rather complicated) bound-state problems 
can be handled "easily," using Mathematica and a few important theorems to calculate 
at least the resulting energy levels to a high precision. We didn't touch upon the task of 
determining also the corresponding wave functions. Not surprisingly, one can show || 
that, for matrix sizes large enough, an arbitrary accuracy for the wave functions can be 
achieved too. On the other hand, for small matrices the general conclusion is that one 
should not trust the wave functions at all, even if the energy levels are rather accurate. 

For the moment, let's only produce three-dimensional plots of our trial function for 
k = and k — 5, to get a feeling how it looks like. The plots are scaled with scale=0 . 7. 



MATHEMATICA 

Reduced trial function ^iy m (x), to be plotted as a function ipk,£m{%, y) of x and y only: 

In[63]:= psix [k_, l_,m_,x_,y_] := Sqrt[2"(2 1+3) k!/Gamma[2 1+3+k]] 
Sqrt [x~2+y~2] "1 Exp [-Sqrt [x~2+y~2] ] LaguerreL[k,2 1+2,2 Sqrt [x~2+y~2] ] 
SphericalHarmonicY [1 ,m,theta,phi] 

Plotting the trial function ipk,£m(%, y) for k = £ = m = 0asa three-dimensional plot. 
This function represents by itself a suitable trial function, namely for the ground state. 



In [64] := Plot3D[psix[0,0,0,x,y] ,{x,-4,4},{y,-4,4}, 
TextStyle->{FontSlant->" Italic" ,FontSize->12}] 
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Plotting the trial function ipk,em{x, y) for k = 5 and I = m = as a three-dimensional 
plot. For k 7^ 0, the functions ipk,£m are not trial functions corresponding to particular 
levels of excitation. Adequate trial functions have to be determined by calculating the 
eigenvectors of the energy matrix; they will be superpositions of various trial functions. 

In[65] := Plot3D [psix [5 , , ,x,y] ,{x,-4,4} ,{y,-4,4} , 
TextStyle->{FontSlant-> " Italic" ,FontSize->12}] 
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